.libPaths('~/R/x86_64-pc-linux-gnu-library/4.1/')
suppressPackageStartupMessages({
library(plotly)
library(gaston)
library(dplyr)
})
# load engine's functions
sapply(FUN = source,
X = list.files("../../src",
pattern = ".R$",
full.names = T))
# R options
options(max.print = 20)
genoFile <- '~/Documents/Work/7_ListenField/genomic/gs-engine/data/check-ddb-simulation/phasedGeno_phased.vcf.gz'
crossTableFile <- '~/Documents/Work/7_ListenField/genomic/gs-engine/data/check-ddb-simulation/crossTable.csv'
SNPcoordFile <- '~/Documents/Work/7_ListenField/genomic/gs-engine/data/check-ddb-simulation/linkageMap.csv'
markerEffectsFile <- '~/Documents/Work/7_ListenField/genomic/gs-engine/data/check-ddb-simulation/markersEffects.json'
predictions <- calc_progenyBlupEstimation(
genoFile = genoFile,
crossTableFile = crossTableFile,
SNPcoordFile = SNPcoordFile,
markerEffectsFile = markerEffectsFile,
# outFile = outFile
)
predictions$Cross <- paste0(predictions$ind1, '_X_', predictions$ind2)
set.seed(1234)
markerEffects <- readMarkerEffects(markerEffectsFile)
crossSimOutFile <- crossingSimulation(
genoFile = genoFile,
crossTableFile = crossTableFile,
SNPcoordFile = SNPcoordFile,
nCross = 100)
simulation <- gaston::read.vcf(crossSimOutFile)
## Warning in gaston::read.vcf(crossSimOutFile): NAs introduced by coercion
## Warning in gaston::read.vcf(crossSimOutFile): Some unknown chromosomes id's (try
## to set convert.chr = FALSE)
simulation <- gaston::as.matrix(simulation)
simulation <- simulation[, row.names(markerEffects$SNPeffects)]
simulation <- data.frame(
simBV = as.vector(simulation %*% as.matrix(markerEffects$SNPeffects)) + markerEffects$intercept,
Cross = as.factor(sub('-\\d+$', '', row.names(simulation)))
)
p <- plot_ly()
# add boxplot of simulated BV
p <- p %>% add_boxplot(data = simulation,
x = ~Cross,
y = ~simBV,
name = "Boxplot of Simulated individuals' BV")
# add jittered markers of simulated BV
p <- p %>% add_boxplot(data = simulation,
x = ~Cross,
y = ~simBV,
line = list(color = 'rgba(0,0,0,0)'),
marker = list(color = 'rgba(31, 119, 180, 0.5)'),
fillcolor = 'rgba(0,0,0,0)',
name = "Simulated individuals' BV",
boxpoints = "all",
pointpos = 0,
jitter = 1,
hoverinfo = 'text',
text = apply(simulation, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
}))
# add predicted BV
p <- p %>% add_markers(
data = predictions,
x = ~Cross,
y = ~blup_exp,
marker = list(color = 'rgb(255,0,0)'),
name = "Theoretical BV",
error_y = ~ list(array = 2*sqrt(blup_var),
type = 'data',
color = '#000000'),
hoverinfo = 'text',
text = apply(predictions, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
})
)
p
alpha <- 0.05
beta <- 0.95
pred_vs_sim <- full_join(simulation, predictions, by = "Cross") %>%
filter(Cross != "Coll0659_X_Coll0425") %>%
group_by(Cross) %>%
summarise(obs_mean = mean(simBV),
blup_exp = unique(blup_exp),
obs_var = var(simBV),
blup_var = unique(blup_var),
p.value = t.test(x = simBV, mu = unique(blup_exp))$p.value,
delta = power.t.test(n = 100, power = beta, sd = sqrt(unique(blup_var)), sig.level = alpha)$delta)
pred_vs_sim$adj.p.val <- p.adjust(pred_vs_sim$p.value, method = "bonferroni")
pred_vs_sim
idLine <- data.frame(mean = c(min(c(pred_vs_sim$blup_exp, pred_vs_sim$obs_mean)) - 1,
max(c(pred_vs_sim$blup_exp, pred_vs_sim$obs_mean)) + 1),
var = c(min(c(pred_vs_sim$blup_var, pred_vs_sim$obs_var)) - 1,
max(c(pred_vs_sim$blup_var, pred_vs_sim$obs_var)) + 1))
rmse <- sqrt(mean((pred_vs_sim$blup_exp - pred_vs_sim$obs_mean)^2))
#plot mean
plot_ly(type = "scatter",
mode = "markers",
data = pred_vs_sim,
x = ~blup_exp,
y = ~obs_mean,
name = "Observed mean vs prediction",
hoverinfo = 'text',
text = apply(pred_vs_sim, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
})
) %>%
add_lines(inherit = FALSE,
x = idLine$mean,
y = idLine$mean,
name = "Identity line")
pred_vs_sim <- full_join(simulation, predictions, by = "Cross") %>%
filter(Cross != "Coll0659_X_Coll0425") %>%
group_by(Cross) %>%
summarise(obs_var = var(simBV),
blup_var = unique(blup_var),
p.value = ks.test(simBV, "pnorm", mean = unique(blup_exp), sd=sqrt(unique(blup_var)))$p.value)
pred_vs_sim$adj.p.val <- p.adjust(pred_vs_sim$p.value, method = "bonferroni")
pred_vs_sim
#plot var
plot_ly(type = "scatter",
mode = "markers",
data = pred_vs_sim,
x = ~blup_var,
y = ~obs_var,
name = "Observed variance vs prediction",
hoverinfo = 'text',
text = apply(pred_vs_sim, 1, function(l) {
paste(names(l), ":", l, collapse = "\n")
})
) %>%
add_lines(inherit = FALSE,
x = idLine$var,
y = idLine$var,
name = "Identity line")
## Document generated in:
## Time difference of 27.34104 secs
##
## CPU: AMD Ryzen 5 3600X 6-Core Processor
## Memory total size: 32.79236 GB
##
##
## Session information:
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Pop!_OS 22.04 LTS
##
## Matrix products: default
## BLAS: /usr/local/lib/R/lib/libRblas.so
## LAPACK: /usr/local/lib/R/lib/libRlapack.so
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.10 gaston_1.5.9 RcppParallel_5.1.6 Rcpp_1.0.9
## [5] plotly_4.10.1 ggplot2_3.4.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 xfun_0.36 bslib_0.4.2
## [4] memuse_4.2-2 purrr_1.0.1 splines_4.1.1
## [7] lattice_0.20-44 colorspace_2.0-3 vctrs_0.5.1
## [10] generics_0.1.3 htmltools_0.5.4 viridisLite_0.4.1
## [13] vcfR_1.13.0 mgcv_1.8-36 yaml_2.3.6
## [16] utf8_1.2.2 rlang_1.0.6 jquerylib_0.1.4
## [19] pillar_1.8.1 glue_1.6.2 withr_2.5.0
## [22] DBI_1.1.3 lifecycle_1.0.3 stringr_1.5.0
## [25] munsell_0.5.0 gtable_0.3.1 htmlwidgets_1.6.1
## [28] breedSimulatR_0.3.2 evaluate_0.20 knitr_1.41
## [31] permute_0.9-7 fastmap_1.1.0 crosstalk_1.2.0
## [34] parallel_4.1.1 fansi_1.0.3 pinfsc50_1.2.0
## [37] scales_1.2.1 cachem_1.0.6 vegan_2.6-4
## [40] jsonlite_1.8.4 digest_0.6.31 stringi_1.7.12
## [43] grid_4.1.1 cli_3.6.0 tools_4.1.1
## [46] magrittr_2.0.3 sass_0.4.4 lazyeval_0.2.2
## [49] tibble_3.1.8 cluster_2.1.2 ape_5.6-2
## [52] tidyr_1.2.1 pkgconfig_2.0.3 ellipsis_0.3.2
## [55] Matrix_1.5-3 MASS_7.3-58.1 data.table_1.14.6
## [58] assertthat_0.2.1 rmarkdown_2.19 httr_1.4.4
## [61] rstudioapi_0.14 R6_2.5.1 nlme_3.1-152
## [64] compiler_4.1.1
shiny::HTML('
<!-- Add icon/font library -->
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.13.0/css/all.min.css">
<link rel="stylesheet" href="https://fonts.googleapis.com/css?family=Lato">
<div class="footer" style="font-family: Lato">
<hr />
<p style="text-align: center;"><a href="https://github.com/juliendiot42">Julien Diot</a></p>
<p style="text-align: center;"><span style="color: #808080;"><em>juliendiot@ut-biomet.org</em></span></p>
<!-- Add font awesome icons -->
<p style="text-align: center;">
<a href="https://github.com/juliendiot42" class="fab fa-github"></a>
<a href="https://www.linkedin.com/in/julien-diot-949592107/" class="fab fa-linkedin"></a>
<a href="https://orcid.org/0000-0002-8738-2034" class="fab fa-orcid"></a>
<a href="https://keybase.io/juliendiot" class="fab fa-keybase"></a>
</p>
</div>
')